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--^7  This  report  discusses  the  problem  of  operating  a  multipurpose  reservoir 
through  regulation  of  a  multilevel  outlet  works  for  a  number  of  water  quality 
objectives.  Operation  of  a  reservoir  to  meet  downstream  goals  for  multiple 
water  quality  parameters  often  results  in  conflict.  A  problem  formulation  and 
solution  are  presented  as  an  attempt  to  reaolve  these  conflicts.  The  multi- 
parameter  reservoir  regulation  problem  is  formulated  in  terms  of  a  scalar 
objective  function,  which  indicates  the  relative  value  of  any  specified  — : :=*. 
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operation  strategy,  and  a  linear  constraint  set.  These  constraints  include 
the  hydraulic  characteristics  of  the  outlet  works  and  any  specified  bounds  on 
the  release  concentrations  of  the  water  quality  parameters.  Two  different 
problem  formulations  are  addressed.  The  target-concentration  problem  is 
formulated  to  achieve  specific  downstream  target  concentrations  without  actual 
constraints  on  the  release  concentrations.  The  constrained-concentration 
problem  is  formulated  to  allow  the  specification  of  upper  and  lower  bounds  for 
all  or  some  of  the  water  quality  constituents.  Both  formulations  can  ac¬ 
curately  deal  with  the  hydraulic  complexity  of  a  multilevel  outlet  works. 

The  algorithms  presented  herein  can  be  used  to  regulate  a  reservoir  in 
a  real-time  mode  in  which  the  state  of  the  system  is  known  by  actual  measure¬ 
ments.  The  algorithms  can  also  be  used  with  an  ecosystem  simulation  model 
in  which  the  state  of  the  system  is  predicted. 
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PREFACE 


The  analysis  and  technique  development  reported  herein  were  con¬ 
ducted  under  Task  I1D.2,  Reservoir  Regulation  Techniques  for  Water 
Quality  Management,  of  the  Environmental  and  Water  Quality  Operational 
Studies  (EWQOS).  The  analysis  and  development  were  conducted  during  the 
period  May  1978-August  1980  in  the  Hydraulics  Laboratory  of  the  U.  S. 
Army  Engineer  Waterways  Experiment  Station  (WES)  under  the  general  di¬ 
rection  of  Mr.  H.  B.  Simmons,  Chief  of  the  Hydraulics  Laboratory, 

Mr.  John  L.  Grace,  Jr.,  Chief  of  the  Hydraulics  Structures  Division,  and 
Dr.  Dennis  R.  Smith,  Chief  of  the  Reservoir  Water  Quality  Branch.  The 
analysis  and  development  were  conducted  and  the  report  prepared  by 
Dr.  Aubrey  B.  Poore  and  Mr.  Bruce  Loftis.  Dr.  Jerome  L.  Mahloch  of  the 
Environmental  Laboratory  was  the  Program  Manager  of  EWQOS. 

Commander  and  Director  of  WES  during  the  preparation  and  publica¬ 
tion  of  this  report  was  COL  Tilford  C.  Creel,  CE.  Mr.  Fred  R.  Brown  was 
Technical  Director. 

This  report  should  be  cited  as  follows: 

Poore,  A.  B. ,  and  Loftis,  B.  1983.  "Water  Quality  Optimiza¬ 
tion  Through  Selective  Withdrawal,"  Technical  Report  E-83-9, 

U.  S.  Army  Engineer  Waterways  Experiment  Station,  Vicksburg, 

Miss . 
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WATER  QUALITY  OPTIMIZATION  THROUGH  SELECTIVE  WITHDRAWAL 


PART  I:  INTRODUCTION 

1.  As  a  result  of  increasing  public  awareness  and  recent  State 
and  Federal  legislation,  water  quality  considerations  in  the  operation 
of  water  resources  systems  are  assuming  a  significant  priority.  Multi¬ 
level  outlet  works  are  the  primary  method  for  controlling  the  quality  of 
releases  from  stratified  lakes.  These  structures  allow  the  release  of 
water  from  various  vertical  strata  in  the  lake.  Typically,  lakes  are 
operated  to  enhance  the  water  quality  of  a  downstream  fishery  or  to  main¬ 
tain  preproject  downstream  temperature  characteristics,  but  a  wide  range 
of  operating  strategies  can  be  envisioned  to  satisfy  different  water 
quality  needs. 


Purpose 

2.  The  purpose  of  this  report  is  to  discuss  the  problem  of  opera¬ 
tion  of  a  multilevel  outlet  structure  for  water  quality  purposes  and  to 
present  algorithms  for  identifying  an  optimal  operating  strategy  that 
considers  many  different  water  quality  constituents. 

3.  The  algorithms  presented  are  for  use  within  a  numerical  model 
that  simulates  the  ecosystem  of  a  lake  through  time  but  could  be  applied 
to  operation  of  a  real-time  system. 

Overview 

4.  Operation  of  multipurpose  lakes  to  meet  water  quality  objec¬ 
tives  is  fundamentally  a  matter  of  making  trade-offs  among  such  aspects 
as  ports,  water  quality  parameters,  time  periods,  and  even  project  pur¬ 
poses  and  projects.  Further,  in  a  density-stratified  lake,  concentra¬ 
tions  of  water  quality  constituents  vary  with  elevation.  Thus,  for 
projects  with  multilevel  intake  structures,  trade-offs  among  ports  must 


be  evaluated  to  determine  which  ports  should  be  opened  and  what  the  flow 
rate  should  be  through  each  port  in  order  to  meet  downstream  water 
quality  requirements. 

5.  Another  type  of  trade-off  in  lake  management  exists  among 
water  quality  parameters.  Often,  operation  to  meet  goals  for  several 
water  quality  parameters  results  in  conflict,  an  example  of  which  is 
downstream  water  quality  requirements  for  water  with  both  a  low  tempera¬ 
ture  and  a  high  dissolved  oxygen  content  to  promote  a  cold-water  fishery. 

In  the  summer,  water  in  the  top  of  the  pool  usually  has  a  high  tempera¬ 
ture  and  high  dissolved  oxygen  (DO)  content,  whereas  water  in  the  bottom 
of  the  pool  is  normally  low  in  temperature  and  DO  concentrations.  A  deci¬ 
sion  must  often  then  be  made  to  either  withdraw  surface  water,  thereby  ful¬ 
filling  the  DO  requirement  but  failing  to  meet  the  temperature  require¬ 
ment,  or  to  withdraw  bottom  water,  thereby  achieving  the  temperature  ob¬ 
jective  but  failing  to  meet  the  DO  objective.  A  possible  trade-off  be¬ 
tween  parameters  would  be  to  assign  weights  to  each  of  them  and  operate 

to  meet  both  objectives  as  closely  as  possible  by  mathematically  mini¬ 
mizing  deviations  from  the  objectives  in  accordance  with  the  specified 
weights . 

6.  A  third  type  of  trade-off  involves  balancing  quality  in  time. 

It  is  possible  to  operate  optimally  on  a  day-to-day  basis  with  no  antici¬ 
pation  of  future  conditions.  Such  operation,  which  can  be  referred  to 

as  static-optimal,  will  result  in  deviations  from  the  objective  through¬ 
out  the  simulation  period.  For  example,  a  lake  can  be  operated  static- 
optimally  to  meet  target  temperatures  in  spring  and  summer.  However, 
because  bottom  water  is  required  to  achieve  the  downstream  target  temper¬ 
ature  in  early  and  mid  summer,  the  cold  water  needed  to  meet  temperature 
targets  in  late  summer  and  fall  can  be  depleted  and  large  deviations 
then  occur.  However,  if  future  conditions  are  known  or  can  be  estimated, 
dynamic-optimal  rather  than  static-optimal  operation  can  smooth  out  the 
daily  deviations  to  achieve  an  overall  lower  level  of  deviation.  The 
principle  behind  dynamic-optimal  operation  is  that  a  large  number  of 
small  deviations  can  have  a  more  favorable  impact  on  the  ecosystem  than 
a  smaller  number  of  larger  deviations.  Dynamic-optimal  operation  can 
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result  in  release  temperatures  that  are  slightly  warmer  than  desired  in 
summer,  thereby  saving  cold  water  for  fall  requirements. 

7.  A  fourth  type  of  trade-off  requires  making  decisions  relating 
to  total  flow  released  from  a  lake  in  order  to  balance  water  quality 
benefits  with  benefits  or  demands  from  other  project  purposes  such  as 
flood  control  or  hydropower  production.  During  a  major  flood  it  is  ap¬ 
parent  that  operation  of  a  lake  for  flood  protection  has  priority.  But 
during  normal  hydrologic  events,  it  is  possible  to  consider  trade-offs 
between  water  quality  and  other  project  purposes.  For  example,  a  peak¬ 
ing  power  operation  could  require  so  large  a  release  flow  rate  that  the 
withdrawal  zone  would  extend  to  the  bottom;  anoxic  water  could  then  be 
released  downstream  or  the  bottom  sediments  could  be  disturbed.  There¬ 
fore,  reducing  the  peaking  power  could  improve  the  release  water  quality 

8.  A  final  type  of  trade-off  concerns  the  operation  of  a  system 
of  lakes  and  connecting  streams  to  achieve  water  quality  requirements  at 
different  locations  throughout  the  system. 
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PART  II:  LITERATURE  REVIEW 


9.  Simulation  models  which  use  optimization  methods  to  identify 
reservoir  operation  strategies  for  satisfying  water  quality  objectives 
are  relatively  new  tools  in  the  field  of  water  resources  systems  anal¬ 
ysis.  They  have  been  used  most  often  for  planning  and  rarely  for  real¬ 
time  operation.  Optimization  methods  have  been  used  to  determine  size 
and  locations  of  selective  withdrawal  intakes  (Loftis  and  Fontane  1976) 
and  to  develop  improved  or  simpler  operation  techniques  (Patterson  et  al. 
1977,  Maynord  et  al.  1978).  Beard  and  Willey  (1970)  developed  a  thermal 
simulation  model  that  included  a  heuristic  procedure  to  anticipate  future 
temperature  objectives  in  determining  reservoir  operation  strategies. 
Kaplan  (1974)  combined  a  lake  ecosystem  simulation  model  and  a  nonlinear 
optimization  technique  to  determine  the  best  operation  of  a  selective 
withdrawal  outlet  structure  considering  constraints  of  several  water 
quality  parameters.  A  scalar  water  quality  index  that  conmensurated  and 
prioritized  several  water  quality  objectives  was  used  as  the  objective 
function  for  this  optimization  problem.  Farber  and  Labadie  (1978)  at 
Colorado  State  University  combined  a  state-space  dynamic  programming 
algorithm  with  the  "WESTEX"  Reservoir  Heat  Budget  Model  (Loftis  1979). 
This  combination  provided  a  systematic  procedure  for  determining  release 
temperature  regulation  strategies  that  ahticiipated  future  meteorological 
and  hydrological  conditions.  Dynamic  programming  was  selected  because  it 
could  handle  sequential  decisions  and  system  nonlinearities  conveniently. 
Fontane,  Labadie,  and  Loftis  (1982)  developed  the  technique  of  objective- 
space  dynamic  programming  which  allocated  violations  of  release  tempera¬ 
ture  from  downstream  target  temperature  such  that  an  objective  function 
for  the  entire  stratification  cycle  was  minimized.  The  Hydrologic  En¬ 
gineering  Center  has  developed  the  model  HEC-5Q  (1981)  to  determine  op¬ 
eration  strategies  for  a  system  of  lakes  and  connecting  streams  consider¬ 
ing  water  quality  as  one  of  many  project  purposes.  The  solution  tech¬ 
nique  presented  in  this  report  is  included  in  the  HEC-5Q  model. 
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PART  III:  A  REVIEW  OF  KAPLAN’S  WORK 


10.  Kaplan  (1974)  presented  one  solution  to  the  multiparameter 
regulation  problem.  He  used  a  lake  ecosystem  simulation  model  (developed 
by  Water  Resources  Engineers  1969)  to  determine  the  water  quality  state 
in  each  of  the  various  layers  of  a  stratified  lake.  From  the  system 
state,  Kaplan  extracted  the  states  of  those  layers  that  contained  ports 
and  constructed  a  state  matrix  $  such  that  the  p  column  represented 
the  water  quality  state  at  the  p  port. 

11.  As  the  objective  function  for  the  optimization  problem,  Kaplan 
used  a  scalar  water  quality  index  which  was  a  function  of  the  release 
concentrations  of  the  water  quality  constituents  under  consideration. 

Upper  and  lower  bounds  on  the  release  concentrations  were  problem  con¬ 
straints.  Kaplan  discovered  at  least  two  serious  difficulties  with  the 
optimization  problem. 

12.  The  first  difficulty  was  that  the  number  of  constraints  was 
at  least  two  or  three  times  the  number  of  independent  variables;  thus, 
it  was  not  always  possible  to  satisfy  all  constraints.  Kaplan  devised 
an  elaborate  scheme  involving  three  sets  of  constraints — most  stringent, 
stringent,  and  least  stringent — to  resolve  this  problem  of  feasible  and 
nonfeasible  constraints  and  the  various  trade-offs.  If  the  most 
stringent  set  of  constraints  could  not  be  satisfied,  the  stringent  set 
or  the  least  stringent  set  was  tried  until  the  constraints  could  be  met. 

A  penalty  was  associated  with  the  stringent  and  least  stringent  sets. 
Kaplan  worked  out  a  trade-off  between  maximizing  the  water  quality  index 
and  minimizing  the  penalty  associated  with  various  constraint  levels  and 
then  decided  on  an  appropriate  set  of  port  openings. 

13.  The  second  difficulty  was  the  reported  existence  of  multiple 
maxima;  thus,  to  find  the  optimum  set  of  port  openings,  all  local  maxima 
points  had  to  be  found  and  the  corresponding  functional  values  compared. 

To  find  as  many  local  maxima  as  possible,  Kaplan  used  a  random  number 
generator  to  generate  starting  points  in  the  feasible  search  region  for 
the  optimization  code.  These  multiple  aiaxima  were  then  used  in  a  trade¬ 
off  analysis  of  constraint  level  and  the  value  of  the  water  quality  index. 
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14.  The  optimization  code  used  by  Kaplan  to  solve  the  optimiza¬ 
tion  problem  was  an  all-purpose,  parameter-free,  penalty  function  code 
developed  at  the  University  of  Texas  by  Staha  and  Himmelblau  (1972). 

This  code,  COMET,  was  developed  to  handle  any  set  of  general  algebraic 
constraints  and  makes  no  allowances  for  a  special  structure  such  as 
linearity.  Kaplan  reported  that  the  code  occasionaly  terminated  due  to 
round-off  error  and  had  to  be  restarted  to  reach  a  local  maximum.  This 
termination  was  accompanied  by  an  objective  function  with  a  flat  surface, 
which  suggests  that  there  may  not,  in  fact,  have  been  multiple  maxima. 

15.  Once  the  local  maxima  were  found,  the  optimal  flow  rate 

through  each  port  was  determined;  then  the  port  openings  were  used  for 
the  next  simulation  period.  The  water  quality  index  was  found  for  each 
layer  of  the  lake  down  to  the  thermocline;  the  arithmetic  mean  of  the 
indices  was  computed  arid  used  as  a  representative  water  quality  index 
for  the  lake  The  total  water  quality  index  was  computed  as 
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where  iu^  and  u)^  are  positive  weights  indicating  the  relative  value 
of  satisfying  in-lake  water  quality  requirements  and  downstream  water 
quality  requirements.  By  comparing  the  WQI  and  the  number  and  severity 
of  the  violated  constraints,  Kaplan  was  able  to  decide  on  a  set  of  port 
openings  to  manage  the  water  quality  of  the  lake  and  river. 

16.  Kaplan  noted  that  Staha  and  Himmelblau  compared  the  COMET  al¬ 
gorithm  to  three  nonlinear  programming  codes  for  25  test  problems.  The 
latter  investigators  found  that  with  analytically  supplied  derivatives 
COMET  was  decidedly  more  efficient  than  the  other  codes,  but  it  was  less 
so  for  cases  where  numerical  approximations  to  the  derivatives  were  used. 
The  derivatives  for  Kaplan's  water  quality  problem  were  determined  numer¬ 
ically,  so  it  is  difficult  to  compare  the  efficiency  of  COMET  with  that 
for  each  of  the  other  nonlinear  codes  tested.  Further,  the  nonlinear 
constraints  can  be  transformed  to  linear  constraints  which  are  handled 
far  more  efficiently  by  primal  methods  than  penalty  function  methods. 
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17.  The  problem  of  multiple  maxima  is  always  difficult  to  resolve. 
The  fact  that  the  objective  function  is  reportedly  flat  suggests  that 
the  algorithm  COMET  stopped  due  to  small  changes  in  the  function  values, 
even  though  the  maximizing  points  were  far  from  the  actual  solution. 


PART  IV:  PROBLEM  FORMULATION 


Water  Quality  Indices 

18.  Construction  of  a  water  quality  index  is  a  mathematical  ap¬ 
proach  that  aggregates  information  on  one  or  more  water  quality  parame¬ 
ters  to  produce  a  single  number  which  indicates  the  relative  quality  of 
the  water  under  consideration.  Such  a  scalar  index  is  essential  for  a 
mathematical  optimization  solution  to  the  problem  of  operation  of  a  lake 
for  water  quality  management.  In  order  to  compute  a  water  quality  index 
for  water  with  known  concentrations,  it  is  first  necessary  to  compute 
subindices  for  each  water  quality  constituent.  For  a  concentration  Yc 
of  a  water  quality  constituent  c  ,  there  is  an  associated  subindex  Sc 
that  measures  the  quality  of  the  water  based  only  on  constituent  c  . 
Graphs  of  subindex  value  versus  concentration  for  several  constituents 
as  suggested  by  the  National  Sanitation  Foundation  (NSF)  are  reprinted 
from  Kaplan's  thesis  with  permission  and  are  presented  in  Figure  Al; 
polynomial  approximations  to  the  subindex  functions  are  presented  in 
Table  Al;  both  can  be  found  in  Appendix  A. 

19.  There  are  several  ways  to  combine  the  water  quality  subindices 
into  a  scalar  index.  The  algorithm  presented  in  this  report  used  an  ad¬ 
ditive  NSF  index 


WQI 


c=l 


(2) 


where 

WQI  =  scalar  water  quality  index 
c  =  index  for  constituents 

=  total  number  of  water  quality  constituents 
iu^  =  relative  weighting  for  constituent  c 
S^.  =  water  quality  subindex  value  for  constituent  c 


The  weights  are  restricted  by 


0  <  ujc  <  1  ;  for  all  c 


(3) 


and 

N 

C=1 


(M 


20.  Another  form  of  the  water  quality  index  is  the  multiplicative 
form  as  presented  by  the  Environmental  Protection  Agency  (Ott  1978). 

N 

_JL  “> 

WQI  =  I]  ScC  (5) 

c=l 

where  the  weights  ui^  satisfy  the  same  restrictions  as  they  do  for  the 
additive  water  quality  index.  The  multiplicative  form  is  more  sensitive 
to  a  low  subindex  value  for  a  particular  constituent  and  could  therefore 
have  advantages  fqr  some  applications.  It  is  more  difficult  to  use  the 
multiplicative  index  for  an  optimization  problem,  however,  because  de¬ 
rivatives  of  the  water  quality  index  for  individual  concentrations  are 
quite  cumbersome. 


Objective  Function 

21.  Two  different  types  of  multiparameter  regulation  problems  can 
be  considered.  The  first,  called  the  constrained-concentration  problem, 
allows  an  acceptable  range  of  release  concentrations  for  each  water  qual¬ 
ity  constituent.  The  second,  called  the  target-concentration  problem, 
is  formulated  to  achieve  specific  downstream  target  concentrations  with¬ 
out  actual  constraints  on  the  release  concentrations.  An  objective  func¬ 
tion  must  be  developed  to  solve  both  types  of  problems.  The  objective 
function  presented  herein  for  each  of  these  problems  is  the  additive 
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form  of  the  water  quality  index;  however,  the  methodology  could  easily 

include  any  other  kind  of  objective  function. 

22.  To  determine  the  value  of  the  objective  function  for  either 

type  of  problem  it  is  first  necessary  to  know  the  state  of  the  system. 

A  system  state  matrix  is  defined  to  contain  the  concentrations  of  the 

various  constituents  at  each  port.  This  concentration  matrix  $  has 

elements  representing  the  characteristic  concentration  of  the 

c*"*1  constituent  at  the  p^  port.  If  q  denotes  the  flow  rate  out  of 
th 

the  p  port,  then  the  release  concentration  for  constituent  c  can  be 
determined  as  a  flow-weighted  average  of  the  concentration  at  each  port. 


2 


Is 

P=1 


where 


c  =  index  for  constituents 

R  =  release  concentration  for  constituent  c 
c 

p  =  index  for  ports 

Np  =  number  of  ports 

q^  =  flow  rate  through  port  p 

*  =  concentration  of  constituent  c  at  port  p 

c,p 

Equation  6  is  not  valid  if  the  water  quality  constituent  under  considera¬ 
tion  is  pH.  One  method  of  computing  the  release  pH  Rc  is 


Rc  =  'l0*l 


I  v  c,p 

P=1 


10  N 


Is 

p=l 


For  the  constrained- concentration  problem,  the  reference  concentration 
Yc  froa  which  the  water  quality  subindex  for  constituent  c  can  be  coa- 
puted  is 

Y  =  R  (8) 

c  c 

For  the  target-concentration  problem, 


Y  =  R  -  T 
c  c  c 


where  equals  downstream  target  concentration  for  constituent  c  . 

Equation  9  is  used  for  teaperature  for  either  of  the  problea  formula¬ 
tions.  The  water  quality  subindices  for  each  constituent  are  com¬ 

puted  as  a  polynomial  function  of  the  reference  concentration  Y 


S  =  S  (Y  ) 
c  c'  c 


For  the  constrained-concentration  problem,  suggested  subindex  functions 
are  determined  by  curve-fitting  the  NSF  functions  in  Figure  Al.  Coeffi¬ 
cients  for  these  functions  are  presented  in  Table  Al.  For  the  target- 
concentration  problea,  the  subindex  functions  can  be  specified  as  para¬ 
bolic,  with  the  reference  concentration  at  the  vertices  of  the  parabolas. 
Suggested  coefficients  for  the  target-concentration  subindex  polynomials 
are  presented  in  Table  A2. 

23.  Thus,  the  subindex  functions  can  be  written  as  follows 


sc  *  2  • 


C  =  1 . N 


where 


c  =  index  for  constituents 

S  =  subindex  value  for  constituent  c 
c 

k  =  counter  for  terms  in  polynomial  representation  of  subindex 
functions 

0  =  highest  order  term  in  polynomial  representation  for  con¬ 

stituent  c 
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a.  =  polynomial  coefficient  for  term  k  and  constituent  c 

K  >  C 

Y  =  reference  concentration  for  constituent  c 
c 

N  =  number  of  constituents 
c 

The  objective  function  to  be  maximized  can  be  written  as 


WQi  =  1  s 

c=l 


(2  bis) 


Constraints 

24.  Operation  of  a  lake  for  water  quality  management  is  con¬ 
strained  by  characteristics  of  the  outlet  structure.  The  hydraulic 
structure  under  consideration  is  illustrated  in  Figure  1.  The  signifi¬ 
cant  components  of  the  structure  are  two  selective  withdrawal  wetwells, 
each  with  a  number  of  selective  withdrawal  ports,  and  a  single  floodgate 
for  larger  releases  independent  of  the  selective  withdrawal  system.  Con¬ 
straints  on  the  outlet  system  include  minimum  and  maximum  flow  rates 
through  each  of  the  ports.  A  hydraulic  constraint  known  as  thermal 
blockage  requires  that  only  one  port  in  each  wetwell  can  be  open  at  any 
time.  Thus,  using  the  floodgate  and  one  port  from  each  wetwell,  a  maxi¬ 
mum  of  three  ports  can  be  open  at  any  one  time.  The  hydraulic  con¬ 
straints  can  be  expressed  as  follows 

F  .  <  q  <  F  5  P  =  1 . N  (12) 

min,p  -  t>  -  max,p  r  p 


where 

F  .  =  minimum  flow  rate  that  can  be  released  through  an  open 

min.p  .  or 

r  port  p 

p  =  index  for  ports 

q^  =  flow  rate  released  through  port  p 

F  =  maximum  flow  rate  that  can  be  released  through  port  p 

max,p 

Np  =  total  number  of  ports  in  outlet  system 

This  relationship  applies  only  to  open  ports.  A  flow  rate  of  zero, 
which  indicates  a  closed  port,  is  also  feasible. 
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Figure  1.  Example  selective  withdrawal  structure 

25.  The  total  flow  through  all  of  the  ports  can  be  constrained  to 
be  (a)  equal  to  a  specified  flow  or  (b)  within  a  range  of  flows.  Either 
constraint  can  be  expressed  as 


where  Qjower  and  Qupper  are  t*ie  minimum  and  maximum  acceptable  total 
flow.  If  Qlowgr  equals  QUpper  >  then  the  flow  constraint  is  an  equality 


constraint. 


26.  The  constrained-concentration  problem  formulation  has  upper 
and  lower  bounds  on  the  release  concentrations  of  each  constituent. 
These  concentration  constraints  take  the  form 


Y.  <  Y  <  Y  ;  c  =  1.....N 

lower, c  -  c  -  upper, c  ’  c 


where 


Y.  =  lower  bound  for  reference  concentration  for  con- 

iower ,c  ... 

stituent  c 

c  =  index  for  constituents 

Y  =  reference  concentration  for  constituent  c 
c 

Y  =  upper  bound  for  reference  concentration  for  con- 

upper,  c  **. 

stituent  c 

Nc  =  number  of  constituents  under  consideration 


Formulation 

27.  The  optimization  problem  can  be  written  to  maximize  a  water 
quality  index  subject  to  hydraulic,  flow,  and  concentration  constraints. 
The  concentration  constraints  are  present  for  the  constrained- 
concentration  problem  and  not  present  for  the  target-concentration  prob¬ 
lem;  the  decision  variables  are  (a)  which  ports  should  be  open  and 
(b)  what  flow  rate  should  pass  through  each  open  port.  Because  there  is 
the  hydraulic  constraint  of  only  one  open  port  per  wetwell,  a  sequence 
of  optimization  problems  must  be  solved.  Each  of  these  problems  has  the 
following  form: 


Maximize 


2  MI 


a.  Yk 
k,c  c 
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subject  to 


nun,p 


-  S  ~ 


roax,p 


Q 


lower  - 


<  Q 

-  ^upper 


f  <  Y  <  Y 

lower, c  -  c  -  upper, c 


=  l,.-.,Nc 


It  should  be  noted  that  the  last  constraint,  given  in  Equation 
be  deleted  for  the  target-concentration  problem. 


(12  bis) 

(13  bis) 

(14  bis) 

14,  should 


PART  V:  SOLUTION  PROCEDURE 


28.  This  part  describes  the  solution  procedure  developed  to  deter¬ 
mine  optimum  strategy  for  multiparameter  reservoir  regulation.  It  pre¬ 
sents  (a)  the  algorithm  formulated  to  address  the  selective  withdrawal 
system  constraints  discussed  in  paragraphs  24-27  and  (b)  the  solution 
technique  for  addressing  the  remainder  of  the  multiparameter  constraints. 

Selective  Withdrawal  Outer  Loop 

29.  One  of  the  constraints  of  the  selective  withdrawal  system  is 
that  only  one  port  in  each  wetwell  can  be  open  at  any  time.  This  con¬ 
straint  can  be  expressed  mathematically  by  introducing  variables  that 
can  only  take  on  the  values  0  or  1.  Techniques  exist  for  solving  prob¬ 
lems  with  0-1  variables,  but  only  at  a  large  cost  in  computer  resources. 
Because  the  problem  size  is  small,  then,  the  most  efficient  procedure 
for  expressing  this  constraint  is  a  simple  enumeration  of  alternatives. 

30.  The  algorithm  proceeds  by  considering  a  sequence  of  problems, 
each  representing  a  different  combination  of  open  ports.  For  each  com¬ 
bination  the  optimal  allocation  of  total  flow  to  ports  is  determined  and 
the  value  of  the  objective  function  is  computed  for  the  optimal  allocation 
of  flows.  The  combination  of  open  ports  with  the  largest  objective  func¬ 
tion  or  water  quality  index  and  its  associated  allocation  of  flows  de¬ 
fines  the  optimal  operation  strategy  for  the  time  period  of  interest. 

31.  There  are  four  different  types  of  combinations  of  open  ports. 
For  one-port  problems,  all  of  the  flow  is  taken  from  a  single  port  and 
the  objective  function  is  computed.  For  two-port  problems,  either  com¬ 
binations  of  one  port  in  each  wetwell  or  combinations  of  a  single  port 
with  the  floodgate  are  considered.  For  three-port  problems,  combina¬ 
tions  of  one  port  in  each  wetwell  and  the  floodgate  are  considered.  The 
total  flow  to  be  released  downstream  is  specified;  if  a  range  of  ac¬ 
ceptable  flows  is  specified,  the  flow  is  treated  as  an  additional  de¬ 
cision  variable  and  the  flow  for  which  the  objective  function  is  maxi¬ 
mized  is  also  determined.  It  should  be  noted  that  if  the  minimum 
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allowable  flow  rate  through  any  port  F  .  is  zero,  then  the  one-port 

min,p 

problem  for  that  port  does  not  have  to  be  included  in  the  enumeration. 

In  fact,  if  all  of  the  minimum  flow  rates  are  zero  only  the  three-port 
combinations  need  to  be  solved. 


Overview  of  Solution  Formulation 


32.  The  discussion  of  the  solution  techniques  and  the  algorithm 
for  solving  the  three-port  combination  flow  allocation  problem  begins  by 
noting  that  maximizing  the  water  quality  index  WQI  is  equivalent  to 
minimizing  the  negative  of  WQI  ,  such  that  f(q)  equals  -WQI(q)  . 

The  solution  techniques  are  then  more  easily  explained  in  terms  of  the 
general  problem 

MINIMIZE  f(q)  (16) 

SUBJECT  TO:  Eq  =  e  (17) 

Aq  >  b  (18) 

where  q  is  an  N^- dimensional  vector  representing  the  flow  rates 


and  Np  represents  the  number  of  ports.  E  is  a 


matrix  of  equality  constraint  coefficients.  For  this  problem  there  is  at 
most  one  equality  constraint,  so  E  is  a  row  vector  and  e  is  a  scalar. 
A  is  a  J  x  Np  matrix  of  inequality  constraint  coefficients,  with  J 
representing  the  number  of  inequality  constraints.  The  inequality  con¬ 
straint  on  the  right-hand  side  (b)  is  a  J-dimensional  column  vector. 

The  rows  in  E  are  assumed  to  be  linearly  independent,  while  those  in 
A  are  not.  In  fact,  the  number  of  rows  in  A  can  be  three  to  five 
times  larger  than  the  size  of  q  ,  so  that  the  rows  of  A  are  necessar¬ 
ily  linearly  dependent.  Since  the  objective  function  f(q)  is  nonlin¬ 
ear,  the  problem  is  a  linearly  constrained  nonlinear  programming  problem. 

33.  The  field  on  nonlinear  programming  is  generally  decisive  in 
choosing  methods  for  both  unconstrained  and  linearly  constrained  prob¬ 
lems.  For  linearly  constrained  nonlinear  programming  problems,  feasible 
direction  methods  (Avriel  1976)  are  generally  accepted  as  the  best 


optimization  techniques.  Within  this  class  of  techniques  the  best 
algorithms  are  the  generalized  gradient  projection  algorithms. 


Input  Requirements 

34.  Data  requiresients  for  the  algorithm  include  information  about 
the  outlet  structure,  the  water  quality  constituents,  and  the  state  of 
the  system.  Hydraulic  data  include  the  number  of  ports  ,  the  mini¬ 
mum  and  maximum  flow  rate  through  each  open  port  F  .  and  F  , 

the  height  above  the  bottom  of  each  port  center  line  ,  and  the  wet- 
well  identifier  of  each  port.  Selective  withdrawal  ports  are  speci¬ 
fied  to  be  in  wetwell  1  or  wetwell  2.  The  floodgate  is  defined  to  be  in 
wetwell  0.  Constituent  information  includes  a  two-character  label  for 
each  constituent  and  the  relative  weights  u»c  for  adding  the  subindices 
in  the  objective  function.  For  the  constrained-concentration  problem 
formulation,  the  upper  and  lower  acceptable  release  bounds  for  each  con¬ 
stituent  and  a  target  temperature  must  be  specified.  For  the  target- 
concentration  problem  formulation,  downstream  target  concentrations  must 
be  specified  for  each  constituent.  The  system  state  is  defined  by  the 
depth  of  the  pool,  the  flow  rate  Q  to  be  released  downstream  or  the 

upper  and  lower  bounds  Q  and  Q,  on  the  downstream  flow  rate, 

r  xupper  lower 

and  the  concentration  matrix  $  containing  the  concentration  of  each 
constituent  at  every  port. 


Initialization  of  Procedure 


35.  To  initialize  the  procedure,  feasible  set  of  flow  rates  q 
must  be  determined  such  that  the  constraints 


Eq°  =  e 


(17  bis) 


Aq°  >  b 

are  satisfied.  As  an  initial  estimate,  q°  is  taken  to  be 
o  /  o  o  o  \T  . 

a  =  »  a2  »  •  t*H  J  where 


(18  bis) 


* 


I 


(Q,  +  Q  )  (F  .  +  F  ) 

o  lower  upper  min,p  max.p 


;  p=l,...,Np  (19) 


2  >  (F  .  +  F  ) 

/ .  min,p  max,p 


This  particular  choice  has  worked  well  in  the  solution  of  the  problem; 
however,  it  does  not  guarantee  that  all  of  the  constraints  Aq°  >  b  are 
satisfied.  If  Aq°  >  b  ,  the  iteration  procedure  begins.  Otherwise,  a 
phase-one  linear  programming  procedure  can  be  used  to  obtain  a  feasible 
point  when  one  exists.  If  no  feasible  region  exists,  the  q°  defined 
in  Equation  19  is  used.  The  explanation  of  how  particular  violated  con¬ 
straints  are  handled  will  be  given  in  paragraphs  41-45  on  projection 
matrices.  The  procedure  starts  with  q°  as  defined  in  Equation  19  be¬ 
cause  most  often  Aq°  >  b  .  The  iteration  then  begins  in  the  interior 
of  the  feasible  region  and  not  on  the  boundary  which  results  from  the 
use  of  the  phase-one  linear  programming  code;  thus,  movement  along  the 
boundary,  which  slows  the  convergence  drastically,  is  avoided  as  often 
as  possible. 

General  Description  of  Solution 

36.  The  general  scheme  of  the  procedure  is  to  iteratively  gener¬ 
ate  a  finite  sequence 

/  k\  k  =  N 

(a  )  (20) 

v  '  k  =  0 

|( 

where  <j  is  accepted  as  the  optimal  set  of  flows  for  the  given  con¬ 
straints  when  one  of  the  following  is  true: 

q  is  a  Karush-Kuhn-Tucker  point 


k  k-1  ,  „ 

a  -  a  <  e! 


f(ak)  -  f(ak_1)  <  £2 


k  =  N 


- 


.i 


where  N  is  some  preassigned  maximum  number  of  iterations  and  and 

are  preassigned  small  numbers.  However,  if  none  of  these  criteria 

2  k+1 

are  satisfied,  an  updated  solution  vector  of  flow  rates  3  is  ob- 

k  k 

tained  from  £  by  choosing  a  search  direction  d  and  then  performing 

K  k 

a  line  search  for  the  minimum  in  the  direction  d  from  q  .  There 
are  several  possible  choices  for  the  search  direction;  e.g.,  Newton's 
direction,  the  negative  gradient  direction,  or  the  more  general  variable 
metric  direction  (Avriel  1976).  Newton's  direction  at  q  is 

^ne  =  -H^(akWak>  (24) 

where  the  gradient  Vf(q)  is  that  column  vector  whose  element  in  the 

itfl  position  is  3f(q)/3q^  and  q^  is  the  i1"*1  element  in  the  vector 

q  .  The  Hessian  matrix  of  f  ,  Hf(q)  ,  is  that  matrix  whose  element  in 
th  th  2  ^ 

the  1  row  and  j  column  is  3  1 (q)/3q^3q^  •  The  negative  gradient  di¬ 
rection  at  q^  is 

-NG  =  'Vf(9k)  (25) 

The  procedure  presented  herein  uses  a  combination  of  Newton's  direction 

R 

and  the  negative  gradient  direction.  Assuming  that  q  satisfies  the 
constraints,  a  projection  matrix  P  is  constructed  to  satisfy 

EPdk  =  0  (26) 

so  that 

E(qk  +  APdk)  =  e  (27) 

A(qk  +  APdk)  >  b  (28) 

for  all  A  in  some  finite  range  [0  ,  A  1  . 

,  Hifl  X 

37.  Initially  dK  is  chosen  to  be  Newton's  direction,  which  is 
to  be  preferred  near  a  solution  in  the  interior  of  the  feasible  region. 

If  dk  and  Pdk  are  nonzero  descent  directions,  then  the  minimum  of 
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k  k 

f(q  +  APd  )  with  respect  to  A  ovpr  |0  ,  A  )  is  located.  The 

max 

corresponding  value  of  A  at  which  this  minimum  occurs  is  denoted  by 

k  k+ 1 

A  ,  and  the  next  iterate  q  is  defined  by 


k+1 

a 


+  AkPdk 


(29) 


If  Newton's  direction  or  the  projected  Newton's  direction  is  not  a  non- 

zero  descent  direction,  d  is  taken  to  be  the  negative  gradient  of  f 

k  k  k  k  k 

at  q  such  that  d  =  -Vf(q  )  .  In  the  event  d  is  zero,  q  is  ac¬ 
cepted  as  the  optimal  flow  solution.  Otherwise,  the  projection  matrix  P 

k  k 

is  again  constructed,  and  f(q  +  APd  )  is  minimized  over  (0  ,  A  ) 
k  k+1  max 

when  Pd  is  nonzero.  Then  q  is  defined  as  above  in  Equation  29. 

The  case 

Pdk  =  p[-Vf(qk)J  =  0  (30) 

k 

is  special  in  that  q  may  be  a  Karush-Kuhn-Tucker  point  and  therefore 

accepted  as  the  optimal  solution.  (See  paragraphs  53-56.)  However,  if 

k  k 

q  is  not  a  Karush-Kuhn-Tucker  point,  P  can  be  modified  so  that  Pd 

J 

is  nonzero  and  q  can  be  defined  as  before. 

38.  The  minimization  of  f(qk  +  APdk)  over  (0  ,  A  ]  is 
called  a  line  search,  and  Davidon's  cubic  interpolatory  scheme  (Walsh 
1975)  is  used  in  this  algorithm.  In  the  sections  to  follow,  these  basic 
ideas  are  expanded  upon  and  various  facets  are  explained. 


Search  Directions  and  Projection  Matrices 


39.  The  point  q  is  assumed  to  satisfy  the  constraints.  The 
unconstrained  search  direction  from  the  point  q  will  be  either  the 
negative  gradient  direction,  Newton's  direction,  or,  as  mentioned  in 
paragraph  36,  one  of  the  variable  metric  directions  (Avriel  1976). 
Newton's  direction  d^g  ,  as  defined  in  Equation  24,  is  actually  ob¬ 
tained  by  solving  the  system  of  equations 

Hf(ak)dNE  =  -Vf(qk)  (31) 
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by  using  Gaussian  elimination  with  scaled  partial  pivoting. 

k  k 

40.  If  the  point  g  satisfies  Ag  >  b  and  there  are  no  equal¬ 
ity  constraints,  then  the  line  search  can  be  performed.  However,  it  is 

often  the  case  that  one  of  the  constraints  (i.e.,  the  r  n  constraint)  is 
k  th 

active  so  that  A  q  -  b  ,  where  A  is  the  r  row  in  A  .  This 
r3  r  r 

may  be  due  to  the  phase-one  linear  programming  technique  used  to  find  a 
starting  point  30  ,  the  existence  of  solution  on  the  boundary,  or  the 
placement  of  a  on  the  boundary  for  the  line  search.  Now  let  R  be 
an  indexing  set  defined  by 

R  =  |r:  Arak  =  br|  (32) 

k 

so  that  R  enumerates  the  set  of  active  constraints.  If  d  is  the 

k 

search  direction,  it  will  often  happen  that  Afd  is  negative  for  some 
reR  .  Then 


Ar(qk  +  Adk)  <  bf  (33) 

for  positive  A  ,  resulting  in  a  violated  constraint.  The  objective 

then  is  to  construct  a  projection  mat',  ix  P  ,  so  that  Pd  is  a  descent 

Ik  ~  k 

direction  producing  a  negative  Vf(o  )  •  Pd  .  Then  A  Pd  >0  when- 

k  k  r 

ever  A  q  =  b  and  EPd  =  0  .  For  A  in  some  finite  interval 
r3  r 

[0  ,  A  J  ,  the  constraints 
max 

E(ak  +  APdk)  =  e  (27  bis) 

and 

A(ak  +  APdk)  >  b  (28  bis) 


will  then  be  satisfied. 

41.  The  projection  matrix  used  here  is  constructed  as 


P 


(34) 
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where  Ijj  is  the  N  x  N  dimensional  identity  matrix,  M  is  an 

P  f 

x  Mp  matrix  (Nm  <  N^)  whose  rows  are  linearly  independent,  and  H 

denotes  the  transpose  of  M  .  The  rows  of  M  are  composed  of  the  lin¬ 
early  independent  rows  of  E  and,  directly  or  indirectly,  those  rows  of 
A(Af)  for  which  Ar<j  =  br  .  (The  explicit  algorithm  for  constructing 
M  is  given  in  paragraph  48.)  An  important  property  of  the  projection 

matrix  P  is  that  A  Pd  =  0  for  any  vector  d  whenever  A  is  a  row 

r  -  -  r 

in  the  matrix  M  or  linearly  dependent  on  rows  in  M  .  Other  important 
properties  of  this  projection  matrix  P  are  that  PP  =  P 


T 

P  =  P 


and  d  =  Pd  +  -  P^d  .  Pd  is  also  orthogonal  to  ^1^  -  P^d 

42.  An  optimal  choice  for  P  is  one  that  would  satisfy 


A  Pd  =  0  if  A  d  <  0 
r  -  r - 


(35) 


and 


A  Pd  >  0  if  A  d  >  0 
r -  r - 


(36) 


That  this  choice  is  not  always  possible  can  be  seen  from  the  relation 

ArPd  =  Ard  Ajl„  -  P\d  (37) 


which  may  be  negative  because  A 


<(\  ■  ?y 


can  be  larger  than  Afd 


In  his  original  paper,  Rosen  (1960)  avoided  this  problem  by  including  in 

the  matrix  M  ,  either  directly  or  indirectly,  all  rows  A  in  A  for 

r 

which  A  q  =  b  .  Of  course,  the  rows  of  E  must  be  included  in  M 
r-1  r  ’ 

so  that  EP  =  0  . 

43.  The  approach  taken  herein  is  to  begin  building  the  matrix  M 
from  the  largest  linearly  independent  set  of  vectors  in  the  collection 
consisting  of  all  rows  in  the  matrix  E  and  all  rows  Af  of  the  matrix 
A  where 


reR  =  jr:  Argk  =  br| 


(38) 
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- 


and 


Ardk  <0  (39) 

Jr  jr 

It  may  happen  that  reR  and  Ad  >  0  ,  but  A  Pd  <  0  .  This  situ- 

ation  can  be  avoided  by  checking  the  sign  of  A  Pd  for  each  reR  .  If 

k  ^ 

A  Pd  <0  for  some  r  ,  the  row  A  is  added  to  the  matrix  M  and  P 
r  ~  r  k 

is  updated.  This  updating  procedure  is  repeated  until  A^Pd  0  for 
all  reR  . 

44.  There  are  several  reasons  why  this  procedure  for  constructing 
P  should  be  efficient  for  this  problem.  First,  it  allows  the  search 
direction  to  point  into  the  region  from  a  boundary  point.  Many  optimiza¬ 
tion  procedures  lose  their  efficiency  for  constrained  problems  because 
they  artificially  require  themselves  to  stay  on  a  boundary  until  the 
Karush-Kuhn-Tucker  conditions  are  checked.  The  small  number  of  inde¬ 
pendent  variables  (flow  rates)  in  this  problem  keeps  the  updating  process 
from  being  repeated  often.  In  fact,  numerical  experimentation  indicates 
that  the  updating  of  P  rarely  occurs  more  than  once,  and  most  often 
not  at  all. 

45.  There  are  two  problems  that  arise  at  this  point.  If  dk  is 

k 

Newton's  direction,  Pd  may  not  be  a  descent  direction.  In  this  case, 

fc  k 

d  is  redefined  to  be  the  negative  gradient  directions  -Vf(q  )  ,  and 

P  is  recalculated  for  this  new  dk  .  If  dk  =  -Vf(qk)  ,  then  Pdk  is 

either  zero  or  a  descent  direction.  That  Pd  is  a  descent  direction 

follows  from  the  identity 

[J-Vf(qk)jTp[J-7f(ak)j  =  [-Vf(ak)]Vp[-Vf(ik)] 

=  £-Pf(qk)J  £-PVf(qk)J  >  0  (40) 

If  Pdk  is  a  nonzero  descent  direction,  a  line  search  is  performed. 

Otherwise  the  Karush-Kuhn-Tucker  conditions  are  checked  for  the  point 
It  k 

q  Either  q  is  a  Karush-Kuhn-Tucker  point,  or  a  row  in  M  can  be 
removed  and  P  updated  so  that  p£-Vf(qk)J  is  a  nonzero  descent 
direction. 
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Violated  Constraints 

46.  At  times  there  are  water  quality  constraints  that  cannot  be 
satisfied  given  the  constraints  on  the  individual  ports  and  total  flow 
rates.  Thus  the  constraint  will  be  violated  at  each  stage  in  the  itera¬ 
tion.  This  is  handled  by  including  in  the  indexing  set  R  those  rows 
k  k  k 

r  for  which  Arq  <  bf  .  Then  Af(q  -  \Pd  )  -  b  will  only  increase, 
with  the  result  that  the  violated  constraint  will  either  improve  or  stay 
the  same,  but  will  not  worsen. 

Construction  of  the  Projection  Matrix 


47.  The  projection  matrix  P  is  constructed  in  stages  from  rows 

k 

in  the  matrices  £  and  A  which  arise  from  the  constraints  Aq  >  b 

and  Eq  =  e  .  The  projection  matrix  is  constructed  only  in  the  event 

(  k 

there  are  equality  constraints  or  when  the  indexing  set  R  =  <r:  A^q 
=  is  not  empty.  The  projection  matrix  P  is  always  defined  as 


p  =  iN  -  mt(mmt)  M 


(34  bis) 


48.  If  there  are  equality  constraints,  M  is  initially  set  by 
M  =  E  which  is  an  £  x  matrix  with  £  linearly  independent  rows. 

If  there  are  no  equality  constraints,  M  is  initially  defined  to  be 

M  =  Af  ,  a  row  vector,  where  r  is  any  index  in  R  for  which  Af  is 
nonzero. 

49.  Several  indexing  sets  are  needed  in  the  construction  of  P  . 
Suppose  the  unconstrained  search  direction  d  and  projection  matrix  P 
are  as  previously  defined.  Then  let 


R2  =  jreR:  ArPdk  <  oj 


RL  =  {reR2:  Af  is  linearly  dependent  on  rows  in  M} 


R3  =  -JreR:  AfPdk  >  oj 


a.  .  VfZ 
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50.  During  the  construction  of  P  ,  the  indexing  sets  R2  ,  RL  , 

and  R3  change;  however,  R  does  not.  The  first  change  in  P  is  aade 
according  to  the  following  algorithm.  Let  P  be  given  by  Equation  34 

where  M  is  either  Af  or  E  .  For  each  rcR2  ,  proceed  according  to 

the  following  steps: 

a.  Calculate  A  P  . 

-  r 

b.  If  A  P  =  0  ,  add  r  to  RL  . 

r 

c .  If  A  P  #  0  ,  add  A  to  H  ,  r  to  RL  ,  and  update  P  . 

51.  After  this  algorithm  has  been  completed,  the  matrix  P  has 

generally  changed  from  the  original  matrix  P  if  R  is  nonempty.  Thus 

there  may  be  an  reR2  for  which  A  Pd*  >  0  for  the  initial  P  ,  but 

kr  ^ 

for  which  A  Pd  <0  for  the  final  P  constructed  as  in  the  above  al- 
r  - 

gorithm.  This  problem  is  rectified  in  the  following  manner.  Given  the 
projection  matrix  P  from  the  above  algorithm,  the  following  is  repeated 

until  A  Pd*  >  0  for  each  reR3  . 

r  *  k 

a.  Calculate  A  Pd  for  each  reR3  . 

k  r  ” 

b.  If  A  Pd*  >  0  ,  leave  rtR3  . 

r  ~k 

c.  If  A  Pd*  <  0  ,  delete  r  from  R3  and  add  r  to  R2  . 

r  -  -  k 

(1)  If  ArPd*  =  0  ,  add  r  to  RL  . 

(2)  If  ArPd*  <  0  ,  add  Af  to  M  and  update  P  . 

52.  At  the  conclusion  of  the  second  algorithm,  the  projection 

k 

matrix  P  has  the  properties  that  A  P  =  0  for  each  rtR2  ,  A  Pd  >  0 

k  r  r 

for  reR3  ,  and  A^Pd  =  0  for  reRL  .  Numerical  testing  shows  that 

the  second  algorithm  is  rarely  executed. 

The  Karush-Kuhn-Tucker  Conditions 

53.  The  Karush-Kuhn-Tucker  conditions  are  checked  whenever 
p£-Vf(q*)j  =  0  .  In  this  case  the  search  direction  is  zero,  and  either 
P  can  be  modified  to  give  a  nonzero  search  direction,  or  <^*  satisfies 
the  Karush-Kuhn-Tucker  conditions  which  are  the  necessary  conditions  for 
a  minimum  of  f  . 

54.  Let 

W  =  (MMT)  M^-Vf(q*)J  (41) 
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where  M  is  an  N  x  N  matrix  (N  <  N  )  and  the  £  x  N  dimensional 
m  p  m  p  p 

E  is  assumed  to  occupy  the  first  £  rows  of  M  .  The  dimension  of  the 
column  vector  W  is  m  .  (The  matrix  M  is  the  same  one  used  in  the 
definition  of  the  projection  matrix  P  in  Equations  34  and  40.)  If  the 
last  m  -  £  entries  in  W  are  less  than  or  equal  to  zero  and 
P^-Vf (qk)J  =  0  ,  then  the  Karush-Kuhn-Tucker  conditions  are  satisfied 
and  the  procedure  stops  with  qk  as  the  optimal  solution. 

55.  If  p[-Vf(qk)]  =  0  and  one  of  the  last  m  -  £  components  of 
W  is  positive,  one  of  the  rows  of  M  is  removed.  Specifically,  if 

VL  :  max  =  r  >  £  +  1  and  >  o\  (42) 

then  the  j1*1  row  of  M  is  deleted  and  the  projection  matrix  P  is  up¬ 
dated  with  this  new  M  .  The  theory  of  Rosen  (1960)  guarantees  then 
that  p£-Vf(qk)]  is  a  nonzero  descent  direction  and  MjP£-Vf(qk)J  is 
positive.  Numerical  experimentation  indicates  that  this  deletion  of  a 
row  from  M  rarely  occurs  because  of  the  present  construction  of  M  , 
thus  negating  the  necessity  of  modifying  P  in  this  manner. 

56.  Because  P  has  been  changed,  there  is  the  possibility  that 
AfP  jj-Vf  (qk)J  <  0  for  some  reRL  whereupon  there  is  a  degeneracy  and 
the  routine  stops.  However,  it  is  almost  always  the  case  that 

A^P  j^-Vf (qk)J  >  0  for  reRL  and  the  1  ine  search  can  be  initiated. 

The  Line  Search  and  Stopping  Criteria 

57.  When  the  line  search  is  initiated,  the  nonzero  descent  direc- 

k  k 

tion  is  Pd  ,  where  P  is  a  projection  matrix  and  d  is  either  the 

negative  gradient  or  Newton's  direction.  Furthermore,  A  Pdk  >  0  for 

k  r 

all  r  for  which  A  q  <  b  .  Now  define  a  unit  vector 

r*  -  r 


i 


and  observe  that  P  ,  and  thus  u  ,  have  been  constructed  so  that 


and 


E(qk  -  Auk)  =  0 


Ar(ak  ♦  Auk)  >  0 


(44) 


(45) 


for  all  reR  .  The  reaaining  rows  in  A  aust  satisfy  A  q  >  b  ,  but 

k  r—  r 

a  negative  value  of  A  u  will  restrict  the  values  of  A  >  0  for  which 

A  (q  +  Au  )  >  b  Thus  the  aaxiaua  value  of  X  ,  X  for  which 
r  k  k  r 

Ar(q  +  Au  )  >  br  ,  or  for  which  a  violated  constraint  will  not  worsen, 
is  defined  by 

I.  (~  :  Vk  <  0  •  br  -  Mk  <  0  • 


A  =  A  . 

■ax  am 


(43) 


k  k 

If  there  are  no  indices  for  which  b  -  A  a  <  0  and  A  u  <  0  ,  then 

T  r  *  r~ 

A _  is  set  equal  to  a  very  large  nuaber ,  such  as  10  . 


max 


58.  The  line  search  that  is  used  to  locate  that  value  of  A  for 


which 


f(gk  +  Akuk)  =  min  -|f(qk  +  Auk)  :  0  <  A  <  ^1MX| 


(46) 


is  Davidon's  cubic  interpolatory  scheme  described  in  detail  by  Walsh 
(1975).  Once  Ak  has  been  approxiswted,  qk+1  is  defined  to  be 


t.kk 
♦  A  u 


(47) 


59.  The  updated  flow  array  qk+1  then  replaces  qk  as  the  solu¬ 
tion  estisMte.  Convergence  for  this  procedure,  and  thus  the  solution 
qk  ,  is  established  by  satisfaction  of  the  conditions  presented  in  para¬ 
graph  36. 
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PART  VI :  SUMMARY 


60.  The  purpose  of  this  report  was  to  discuss  the  problem  of 
operating  a  multipurpose  reservoir  through  regulation  of  a  multilevel 
outlet  works  for  a  number  of  water  quality  objectives.  Operation  of  a 
reservoir  to  meet  downstream  goals  for  multiple  water  quality  parameters 
often  produces  conflicts;  a  problem  formulation  and  solution  were  pre¬ 
sented  as  an  attempt  to  resolve  these  conflicts.  The  multiparameter 
reservoir  regulation  problem  was  formulated  in  terms  of  (a)  a  scalar  ob¬ 
jective  function  which  indicates  the  relative  value  of  any  specified 
operation  strategy  and  (b)  a  linear  constraint  set.  These  constraints 
include  the  hydraulic  characteristics  of  the  outlet  works  and  any  speci¬ 
fied  bounds  on  the  release  concentrations  of  the  water  quality  parame¬ 
ters.  Two  different  problem  formulations  were  addressed.  The  target- 
concentration  problem  was  formulated  to  achieve  specific  downstream  tar¬ 
get  concentrations  without  actual  constraints  on  the  release  concentra¬ 
tions.  The  constrained-concentration  problem  was  formulated  to  allow 
the  specification  of  upper  and  lower  bounds  for  all  or  some  of  the  water 
quality  constituents.  Both  formulations  can  accurately  deal  with  the 
hydraulic  complexity  of  a  multilevel  outlet  works. 

61.  The  algorithms  presented  herein  can  be  used  in  a  real-time 
mode  in  which  the  state  of  the  system  is  known  by  real-time  measurements 
The  algorithms  can  also  be  used  with  an  ecosystem  simulation  model  in 
which  the  state  of  the  system  is  predicted. 

62.  The  algorithms  provided  an  efficient  procedure  for  solving 
the  multiparameter  reservoir  regulation  problem.  The  linear  structure 
of  the  constraint  set  has  been  used  to  advantage.  The  problem  as  pre¬ 
sented  was  small;  most  of  the  matrices  were  3X3,  which  was  the  dimen¬ 
sion  of  the  number  of  decision  variables  (open  ports).  The  constraint 
matrix  was  larger,  but  calculations  were  done  with  only  those  con¬ 
straints  that  were  active,  which  was  often  just  one  or  two.  This  solu¬ 
tion  procedure  is  particular  to  the  multiparameter  reservoir  regulation 
problem;  thus,  the  size  and  complexity  of  a  general  purpose  nonlinear 
optimization  code  has  been  avoided. 
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63.  The  concept  of  a  scalar  water  quality  index  was  presented. 

It  must  be  emphasized  that  the  index  was  used  as  an  example  objective 
function  for  the  optimization  problem  under  consideration  and  that  this 
is  not  a  recommendation  for  a  particular  water  quality  index  or  even  for 
the  general  concept  of  a  water  quality  index.  It  was  a  useful  tool  for 
presenting  the  optimization  algorithms  because  (a)  it  is  a  single  number 
with  functional  dependence  on  the  various  parameter  concentrations  being 
considered  and  (b)  the  necessary  derivatives  can  be  determined  analyti¬ 
cally.  For  a  particular  application  of  the  optimization  procedure,  a 
much  simpler  objective  function  might  be  entirely  adequate  and  thus  more 
appropriate . 
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APPENDIX  A:  WATER  QUALITY  SUBINDEX  FUNCTIONS 
AND  COEFFICIENTS  (FROM  KAPLAN  1974) 


Table  A2 

Coefficients  for  Target-Concentration 
Subindex  Polynomials 


2 

Polynomials:  f  =  a  +  bx  +  cx 

Where  x  is  Release  Concentration  Minus  Target  Concentration 


Parameters  x 

a 

b 

C 

Temperature  (T  ) 

100 

0 

-4.0 

Acidity  (pH) 

100 

0 

-11.11 

Dissolved  oxygen  (DO) 

100 

0 

-4.0 

Total  solids  (TDS) 

100 

0 

-0.000625 

Biochemical  oxygen  demand  (BOD^) 

100 

0 

-0.444 

Fecal  coliforms  (FColi) 

100 

0 

-0.0025 

Nitrogen  (NO^) 

100 

0 

-0.16 

Phosphorus  (PO^) 

100 

0 

-16.0 
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In  accordance  with  letter  from  DAEN-RDC,  DAEN-ASI  dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a  facsiaile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 
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